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Simulation results of the thermal conductivity £ of Dissipative Particle Dynamics 
model with Energy Conservation (DPDE) are reported. We also present an analysis of the 
transport equations and the transport coefficients for DPDE based on a local equilibrium 
approximation. This approach is valid when the particle-particle thermal conductivity 
A and the friction coefficient ( are large. A qualitative derivation of the scaling form 
of the kinetic contribution of the transport of energy is derived, yielding two different 
forms for the kinetic contribution to the heat transport, depending on the value of A. 
We find agreement between the theoretically predicted value for £ and the simulation 
results, for large A and many particles interacting at one time. Significant differences are 
found for small number of interacting particles, even with large A. For smaller values of 
A, the obtained macroscopic thermal conductivity is dominated by diffusive transport, in 
agreement with the proposed scaling form. 



I. INTRODUCTION 



The model for the simulation of the dynamics of complex systems, known as Dissipative Particle 
Dynamics (DPD), has undergone an important development in recent years since it was introduced in 
1992 by Hoogerbrugge et al. |Q. This off-lattice simulation methodology is especially addressed at the 
modeling of the behavior of fluids at a scale where fluctuations are important. Thus, the method is suitable 
for both the study of the dynamics of small systems and in situations dominated by Brownian motion, for 
instance. The potential of this methodology was realized from the beginning and, as a consequence, DPD 
has already been applied to the study of the dynamics of systems of practical interest such as polymeric 
materials ^-^], polymer adsorption Rj, colloidal suspensions ||, multicomponent systems ||, among 
others. Despite the advances carried out towards a complete understanding of the features of this new 
methodology, there are still unclear aspects, mainly concerned with the large-scale transport properties 
of the model. The study of these properties, relative to the extended version of the DPD algorithm which 
includes heat transport is the subject of this paper. 

In the DPD model, an ensemble of particles, viewed as mesoscopic representations of flocks of physical 
molecules, interact with each other through conservative as well as dissipative forces. Random forces are 
also explicitly considered to maintain the particles in constant agitation even in the absence of exter- 
nal force fields. In this method, the particle-particle interactions are such that the total momentum is 
conserved despite the presence of dissipative and random forces and, therefore, the macroscopic hydro- 
dynamic behavior at long-wavelength and long-time scales is guaranteed. In addition, since the model is 
not defined on a lattice, it is Galilean invariant and has no extra conservation laws apart from those that 
are physically relevant. Together with the theoretical advantages, the DPD simulation methodology can 
be implemented using only local interactions, which makes it fast to run on a computer, as in Lattice 



Gas cellular automata |10|, for instance 



The original DPD method was only concerned with the conservation of the total momentum of the 
particles so that the model could reproduce continuity and Navier-Stokes equations at the macroscopic 
level. Several refinements, such as to demand that the equilibrium momentum distribution of the particles 
be of the Maxwell-Boltzmann type pi] ], were added aiming at a closer relationship of the model with 
real physical systems. However, the energy transport and, therefore, the heat flow, could not be modeled 
since the DPD model was essentially isothermal. 

An extention of the model to also incorporate the conservation of the total energy in the particle-particle 
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interaction, has been proposed in ref. (2) , and later in ref. JT^] . More recently, we have refined the original 
algorithm to ensure energy conservation at every time-step, instead of in the mean ||. In our version of 
the DPD model with energy conservation (referred to as DPDE from now on), together with the position 
and momentum, an internal energy of each particle is explicitly considered as a relevant variable. With 
the definition of a particle's internal energy, an energy balance can be established so that the work done 
by the dissipative and random forces is stored or released from the particle's internal energy, and thus 
the total energy is held constant in the particle-particle interaction. Therefore, as a consequence, the 
addition of the energy conservation to the momentum conservation in the particle-particle interaction 
allows the DPDE model to describe heat transport together with the momentum transport. Hence, the 
five hydrodynamic fields, density, momenta and energy, issued from the conservation laws at a microscopic 
level, can be properly reproduced by the DPDE model. 

In this paper we address the study of the transport properties of the DPDE model from both, theoretical 
as well as simulation points of view. In the case of the isothermal DPD, much work has been devoted 
to the analysis of the transport properties of the model, especially with regard to the calculation of 
the shear viscosity fl3|-|l5|l, towards a qualitative and quantitative understanding of the long- wavelength 
long-time dynamic behavior of the model. Significant progress has been done in this direction by the 
formulation of Boltzmann-type equations, which have been solved by means of different approximations. 
In particular, in ref. jl3| the transport equations for the macroscopic fields were derived and expressions 
for the transport coefficients were given, together with a lucid discussion of the conceptual implications 
of the approach. The authors distinguish a kinetic and a dissipative contribution to the shear viscosity, 
essentially proportional to 1/7 and 7, respectively, with 7 being the friction coefficient per unit of mass of 
the particle-particle dissipative force, in that reference. Roughly, the first one is due to the transport of 
momentum due to the random motion of the particles while the second is the contribution due to the direct 
frictional force between particles. The qualitative behavior of the viscosity coefficient was captured by the 
approach, but significant deviations between the simulation and the theoretical predictions were found. 
The disagreement was found to be more important in the region dominated by the kinetic contribution, 
although theoretical and simulation results seem to tend to agree in the limit where a given particle 
interacts with many others at one time. Such a discrepancy has been further reported by the simulations 
of Pagonabarraga et al. |l6| . However, very recently, the approximations used in the solution of the 
Boltzmann-like equation of ref. |l3) have been reviewed in ref. flq] . Masters et al. claim that the 
discrepancies between the theory of ref. |l3| and the simulations can be explained if the correlations 
between the collisions are taken into account. This is even more dramatic in the case that only a few 
particles are in simultaneous interaction with a given one. 

For the DPDE model, such a deep analysis of the transport phenomena displayed by the system is still 
lacking. In ref. || we analyzed the equilibrium properties of the model from both the simulation and 
the theoretical perspectives. As far as the transport properties of the model are concerned, in ref. || we 
derived, from simple arguments, expressions for the macroscopic thermal conductivity and the viscosity 
coefficients of the DPDE model, in the limit where the dissipative interactions at the mesoscopic level are 
dominant. Furthermore, we have shown figures with convection rolls and temperature profiles in a DPDE 
system with cross temperature gradient and gravity field, proving qualitatively that the model displays 
the same features as expected for a fluid. However, simulation results on these transport coefficients is 
scarce. Only data from a model with particles at rest exchanging heat to each other are available |p7| . 

In the present paper, we do not attempt to derive a complete theoretical treatment of the transport 
properties of the DPDE model, more complex to handle than the isothermal DPD model due to, first, a 
larger number of parameters, second, the couplings between the different macroscopic fields, and third, 
the dependence of the transport coefficients with respect to the temperature. Instead, we propose a 
calculation of the complete set of the transport equations for the five hydrodynamic fields in the limit 
in which the local equilibrium hypothesis is reasonable. Such a procedure is well known for dissipative 



systems [18 20 and has also been discussed in this context in ref. JL3|. This treatment permits us to 
derive transport equations for the hydrodynamic fields, and approximate expressions for the transport 
coefficients, valid in the region of parameters where dissipative contributions are dominant. We report 
simulation data on the thermal conductivity, characteristic of the DPDE model, for different values of the 
parameter describing the direct heat transfer between particles, and for two values of the particle-particle 
interaction range. In all the cases studied, the system can be considered as non-dilute, in the sense that 
there are always many particles interacting with each other at one time. Different behaviors are observed, 
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varying with both, the magnitude of the dissipative coefficient, A in our model, and the interaction number, 
n that we define as the number of particles interacting at one time. We attempt a qualitative explanation 
of these different behaviors. Finally, we also report simulation data of the viscosity of the DPDE model 
as a function of the interaction range. The agreement between the simulation results and the theoretical 
predictions for the transport coefficients is good in the region dominated by the dissipative interaction, 
and for large interaction number, n. We observe, however, significant deviations if the interaction number 
is decreased, regardless of the value of the dissipative coefficients. The theoretical explanation of this fact 
could be related to the recent analysis of ref. |Q for the isothermal DPD model. In the low A limit, we 
observe a particular behavior that we relate to a form of transport due to particle diffusion. 

The paper is organized as follows. In section 2, we give a brief review of the algorithm describing the 
DPDE model used in the simulations. In section 3, we perform the derivation of the transport equations, 
based on a local equilibrium formulation. In this section we give expressions of the transport coefficients 
and thermodynamic properties that are compared with those of ref. ||. In section 4, we show the results 
obtained for the thermal conductivity from the simulations performed with the DPDE model in different 
conditions, and compare with the theoretical predictions. Finally, in section 5 we review and discuss the 
main results of this paper. 



II. THE DPDE MODEL 



The DPDE model is defined from the stochastic differential equations for the updating of the relevant 
variables of the particles |HJ!|. This implies equations for the change in positions, {r"i}, and momenta, {pi}, 
of all the particles, according to Newton's equations of motion. These interactions include frictional and 
random forces due to the mesoscopic nature of the particles, issued from a coarse-grained picture of the 
physical system represented by the DPDE model. In addition, the internal energies {ui} of the particles 
are also considered as relevant variables. The particle internal energy is necessary in this description 
to impose the conservation of the total energy in the particle-particle interaction, when dissipative and 
random forces exist. From a physical point of view, both the random interactions and the frictional 
forces, as well as the internal energy of the particles, arise due to the fact that the DPD particles may 
be thought of as if they were subsystems of a large number of degrees of freedom that are not explicitly 
taken into consideration in the simulation. 

In ref. |3| we have discussed at length the derivation, the properties and physical consistency of the 
algorithms defining the DPDE model. The reader is addressed to that reference for the details. Thus, let 
us consider the dynamics of the DPDE model as that of an ensemble of N particles interacting through 
conservative, dissipative and random interactions. Let Ff xt be the external force acting on the \ th particle, 
and F@, the particle-particle conservative force between the pair ij. The frictional force exerted by the 
j th particle on the i th one is given by 



which is directed along the vector ry = r^/Vy, with fy = (r) — fi) and ry = \ fj — f\\, so that the total 
angular momentum is also conserved. Cij^ij^ij is the particle-particle friction tensor, which is a function 
of the distance between the particles and also of the particles' temperatures 0i, defined later on, through 
the scalar friction coefficient £y = Qj {r ij ,9i,9j). 

As in the isothermal DPD model, we will also demand that the particles arrive at a given thermal 
equilibrium for the momenta, described by the Maxwell-Boltzmann distribution. Thus, together with the 
dissipative forces, we have to also introduce random Brownian forces to compensate, on average, the loss 
of energy due to the frictional forces. These random forces are Gaussian and white, and their amplitudes 



are obtained from a suitable fluctuation-dissipation theorem [21 



As mentioned, the dissipative particles, as mesoscopic views of physical systems, bear a large number 
of internal degrees of freedom. We assume that these degrees of freedom relax fast towards a thermo- 
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dynamic equilibrium characterized by the particles' internal energy itj. This permits us to introduce a 
thermodynamic description at the particle level and thus define a particle's temperature, 9i, as well as a 
particle's entropy, s{ui), related to the internal energy by a particle's equation of state 

1 dstuA ,„ „. 

— = — — (2.2) 

9i du, K ' 



Dynamically, the i mesoscopic particle varies its internal energy due to the work done by frictional and 
Brownian forces on the other particles. However, we consider that it can also vary its internal energy by 
exchanging mesoscopic heat, q®, with the j th one, if they have different particle temperatures: 

Q?j = A« (Oj - 9i) (2.3) 



where Ay = Ay (r.y, 0j, 0j) is the mesoscopic thermal conductivity. Eq. ( |2.3| ) recalls the macroscopic 
Fourier's Law for heat transport, here defined at a mesoscopic level. Therefore, the change of the internal 
energy of a given particle is due to both the mesoscopic heat exchange as well as to the work done by the 
dissipative forces. In addition, we demand that the internal energy be a fluctuating variable described 
by classical equilibrium statistical mechanics |^| as that of a small system in contact with a heat bath. 
In our model, the dynamics of such fluctuations is introduced through a random heat flow between the 
particles, q^ 



that satisfy a given fluctuation-dissipation theorem |^,[ 



The algorithm for the DPDE model described so far is obtained after integration of the stochastic 
differential equations for the change in the relevant variables (Langevin-like equations) in a small time 
step St, and retaining terms up to first order ||. We have chosen an interpretation rule jl8| for these 
stochastic differential equations that is neither Ito nor Stratonovich, but the resulting stochastic process 
satisfies detailed balance, a property required for the system to approach to the proper thermal equilibrium 
7f|,f3|J3l ■ One arrives at an implicit algorithm given by H 



r! = n + !±6t 
m 



Pi = Pi 




u' 4 = Ui + 



E{^[(w-p*o-^] 2 + A y (^-^ 



)\5t 



(2.4) 
(2.5) 



E {^W • ^ ^ + Sign(* - j)Ki ^} St 1/2 



(2.6) 



where St is the time-step. Here, the prime symbol stands for the calculation of the respective variables at 
an advanced instant of time t + St p3). Thus, f/ = fi(t + St), p/ = p>i(t + St), and u\ = Ui(t + St) while 
Ti, pi and Ui are the value of these functions at time t. T'^ and A^- are the amplitudes of the random 
force and heat flow given by the corresponding fluctuation-dissipation theorems H 

rf i = 2fce iJ c« (2.7) 

A?- = 2kXij6i0j (2.8) 



where we have defined the mean inverse temperature as Gy = (l/0j + l/0j)/2. Again, the prime indicates 
that these amplitudes need to be calculated with their arguments evaluated at the advanced time. In eqs. 
(2.5) and (2.6) we have introduced the random numbers ClW and fij* , which are Gaussian and white, 

with correlations (^foff) = (fi^Oj^) = {S lk Sji + S u S jk ) and (fi^O^f ) = 0. These random numbers 
arise from the contributions due to the random forces and random heat flows, respectively. In ref. J3| an 
explicit algorithm for the DPDE model is also given. 



4 



Although the theoretical predictions are valid for any DPDE model, for the simulations done in the 
present work, however, we have used the model defined by the following functions. For each particle 
we have used the equation of state u, = tfiOi , where <p is a constant that stands for the heat capacity of 
a single particle. Notice that it is convenient to choose 3> k to avoid negative values of the internal 
energy of the particle during the simulations. For the dissipative interactions we have chosen 



, 2 

TV 



Cij = Co ( 1 - for rij < r c ; and if r 2j > r c (2.9) 

x v = Tnr f 1 ~ — ) for r « ^ r ^ and if r a > r * ( 2 - 10 ) 



'A 



where Co and £o are constants giving the magnitude of the mesoscopic friction and thermal conductivity, 
respectively, and and r\ are the respective ranges of these dissipative interactions. Note that Ajj 
explicitly depends on the particles' temperatures. Here, no potential interaction forces are considered. 



Despite the multiple parameters that the algorithm given in eqs. ( |2.4| ), (2.5) and ( f2.q ) displays, we use 
dimensionless variables that reduce the relevant parameters to two 

B= k - (2.11) 

O-T^T (2-12) 
Co^ 



_B measures the relative magnitude of the energy fluctuations with respect to the average energy of 
each particle, while C is a ratio between the characteristic time of momentum fluctuations and energy 
fluctuations for one particle. To make the variables dimensionless we have defined a temperature of 
reference Tr to be used as the scale of temperature. Thus, T = TrT* and 9{ = Tr9*, where the 
asterisk is used to denote a dimensionless variable from now on. The momentum of the particles is made 
dimensionless according with pi = p*\f2mkTn, using the characteristic value of the momentum thermal 
fluctuations \J {pf). The momentum penetration depth for a pair of particles defines a characteristic 
length scale I = y / 2mkTfi/(o . Thus, r*j = r*l. The characteristic scale of time is the relaxation time 
for the particle's momentum, that is t = t*ra/C,o- The characteristic scale for the particle's internal 
energy is chosen to be 4>Tr so that Ui = u*4>Tji. In summary, when these dimensionless variables are 
introduced into the algorithm, the equation for the momentum and the position remarkably contains no 
parameters (in the absence of conservative interactions) , and only depends on the dimensionless range of 
the frictional force rt. In turn, the equation of the energy is such that the terms related to the mechanical 
energy dissipation are proportional to B, while the particle-particle heat transport is proportional to C. 
In addition, the equation is also functionally dependent on the ranges rJ and 



III. MACROSCOPIC PROPERTIES OF THE DPDE 



The model described so far has a hydrodynamic behavior at long-wavelengths and long-times, according 
to the conservation laws introduced in the particle-particle interaction. In this section we will derive the 
macroscopic transport equations of the model in the local equilibrium approximation, in order to relate 
the model parameters with transport coefficients and thermodynamic relations. Such a procedure yields 
the dominant contribution of the transport coefficients in the case in which the dissipative interactions 
are large, so that the kinetic transport is subdominant. Under these conditions, we can expect the system 
to relax fast to the local equilibrium probability distribution function for the set of variables describing 
the state of the system, i.e. 

PUin}, { Pl l R}, t) ~ exp j- £ y - + ___ + -g w J 
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(3.1) 



We have used the shorthand notation Vi = v(fi, t), Tj = T(fi, t) to indicate the macroscopic velocity and 
temperature fields, respectively, at the space points occupied by the i th particle. Furthermore, in the 
term involving the particle-particle interaction potential ip we have written T^ 1 = (1/Tj + l/Tj)/2 to 
maintain the symmetry of the probability distribution under the exchange i — ► j. Corrections to this 
probability distribution are essentially proportional to l/£ and to 1/A, and are neglected in this approach 



Let us consider a single component system and ignore external fields for simplicity. The particle number 
density is described by the field in the phase space 



JY 



(3.2) 



\i=l 



where here the average must be taken by means of the instantaneous non-equilibrium probability distri- 
bution, that we will approximate by eq. (3.1). After differentiating both sides of eq. (3.2) and averaging, 
we obtain a continuity equation reflecting the microscopic particle number conservation 



0_ 

at 



P(r,t) 



I N _ 

V — 

m 



5{r 



-V-v(f,t)p(r,t) 



(3.3) 



where, in deriving the first equality of this last equation, use has been made of eq. (2.4) written in 
differential form. 



The momentum density in the system is given by the expression 

/ n \ 



l{r,t) = ( ^2p t S{r- fi) ) = mp(f,t)v(f,t) 



(3.4) 



where the second equ ality is in fact a definition of the baricentric velocity |Bj| v(f,t). Again, time 
diffe rent iating eq.(3.4) and performing the average with the local-equilibrium probability distribution, 
eq. ( |3.l| ), after some algebra one arrives at 



d_ 

dl 



mpv = —V • mvvp — VkT p 



dfi 



-Af- 



dAr 



+ C(Ar, r, r/)AfAf • (v(f/) — v{f)) 



p (2) (r, rl) 



(3.5) 



where we have omitted the explicit time and space-dependence of the macroscopic fields where no con- 
fusion could occur. Here, Ar = rl — r, Ar = |Ar|, and Af = Ar/Ar. In addition, we have defined the 
two-particle density from the relation j25|] 



p W{r,?i,t) = J2 (S(f-n)S(f/-f 3 )) = p(r,t)p(r/,t)g(Ar,r,r/) 



(3.6) 



where the second equality is, in fact, a definition of the pair distribution function g(Ar,r,rf). For this 
function we have explicitly written its dependence with respect to Ar due to the interparticle potential, 
together with t he f and r7-dependence due to the possible spatial variation of the temperature field 
occuring in eq. (3.1). It is important to realize that g{Ar,r,rl) is invariant under the exchange f — > rl. 



G 



Furthermore, in eq. (3J5) we have introduced the energy-averaged friction coefficient, resulting from a 
possible temperature dependence in the mesoscopic friction £, according to 



C(Ar,r,r/) = — 



duiduj Cy (Ar; 9^), ^M/k-u./kT^sM/k-u./kT, 



(3.7) 



where A is a normalization factor. Since all the pairs of particles are equivalent, we have dropped the 
subscript ij on the left hand side of this last equation. Again, the explicit r and r/-dependences come 
from the spatial variation of the macroscopic temperature field and is also inv aria nt under the exchange 



f — > ri by construction of £. To obtain the long-wavelength behavior of eq. (3.5), we use the fact that 



the integrand in eq. (3.5) is significantly different from zero only in a small neighborhood of ri around 
rl = r, because £ is a function of short range with respect to Ar. We thus expand the ri dependence of 
the velocity, temperature and the density fields in powers of Ar, up to the first significant order. Then, 



after some algebra, eq. (3.5) can be cast under the form 

d _ =! 

—mpv = — V • mvvp — Vp + V • a (3-8) 

which is that of a Navier-Stokes equation. Here, we have defined the pressure exerted by the DPDE 
system 

2 -2 r, A .u^) 



p(f, t) = kTp(r, t) - — / / dAr Ar 6 ^-^^(Ar, r, r) (3.9) 
3 J dAr 



Note that in this expression the last argument of g(Ar, r, r) is r and not ri. Eq. (3.9) is formally identical 
to that obtained from thermodynamic considerations for a physical system with pair interaction potentials 
p5[ . A similar expression has been also obtained in ref . for the isothermal DPD model, and is given 
in eq.(3.10) of ref. JU for the DPDE model. We have also defined the macroscopic stress tensor field 

— 5 =: \ 
2r/ iP 2 Vv+-r/ 4 p 2 lV ■ vj (3.10) 



From the coefficient of the symmetric and traceless velocity gradient, Vu, we identify the shear viscosity 
coefficient 

2n f°° 

r) = p r/4 with 774= — / dAr Ar £(Ar,r,r)g(Ar,r,r) (3-11) 
15 Jo 



The coefficient 774 is a function of the position and of the time through the temperature field. The volume 
viscosity is obtained from the coefficient of the V • w-term, and reads 

vv = ^mp 2 (3.12) 



Both viscosity coefficients agree with the dominant contribution of the transport coefficients given in ref. 
|jl3|| , in the limit of negligible kinetic momentum transport. Here, however, these viscosity coefficients can 
be functions of the temperature through £ due to the possible dependence of £ on the particles' temper- 
atures. Note that neither kinetic nor particle-particle interaction potential contributions are significant 
to the momentum transport coefficients in the large friction limit developed here. 

The calculation of the energy transport equation requires the calculation of the macroscopic equations 
for the transport of each of the three contributions to the total energy density, that is, the kinetic, 
potential, and internal energy densities. Let us first analyze the change in the mechanical energy density, 
given by the sum of the first two mentioned contributions. The calculation follows the same lines as in 
the case of the momentum transport, that is, we first time-differentiate the proper variable, 
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p(f,t)e m (f,t) = ^ 



Pt 

2m ' 2 



5(r - fj) 



(3.13) 



in this case, and then expand the result in the gradients of the macroscopic fields. One finally arrives at 
^pe m = -V ■ (^mpv 2 + + 7}P kT ) v-V -vp + v- • ffj (3.14) 



where 



*(f,t) = 2-irp / dAr Ar 2 g(Ar)ip(Ar) 



(3.15) 



is the macroscopic potential energy density. The behavior described by eq. (3.14) was already found in 
ref. |l^] for the long-time behavior of the mechanical energy of the isothermal DPD model. 

The equation describing the transport of the particle's internal energy is found by analyzing the evo- 
lution of the density 



(3.16) 



Time-differentiating this variable and averaging, one finds 
d 



dt 



pe u = — V • pe u v + / dfl 



' [(fj - pi) ■ Arf C(Ar, 9 h 9j) + A(Ar, 9 i ,9 j ){9 j - 9 t ) 



2m 



^ 6(f-fi)S(f/-fj) 



(3.17) 



As before, the integrand is only different from zero for small values of Ar, which allows us to expand all 
of the r/-dependence of the integrand around f. Up to the first significant order we obtain the transport 
equation for this form of energy 



d i 
— pe u = — V ■ pe u v + 3 : Vv + V ■ £VT 
at 



(3.18) 



The first term on the right hand side of eq. ( 3.1 8| ) is the particle internal energy advected by the mean 
flow, the second term is a source of particle internal energy due to the so-called viscous heating due to 
the mechanical energy dissipated by viscous forces. The last term in this equation is the macroscopic 
energy transport due to heat flow. Let us analyze this heat flow term in more depth. First of all, we 
define the energy-averaged function A, according to the relation 

\(Ar,f,ft)(T(ft)-T(f)) = -L j du.duj A(Ar; 6 t , 9 j )(9 j - 9 t ) e «(«i)/*-«*/*r <e ii(uj)/*-i*j/*ai (3.19) 



as in eq. ({Tt|). Thus, from eq. ( [3.17 ) and eq. ( |3.19 ) we can write 

J dft/\(Ar, Oi, 9j)(9j - 0<) S(f- f) 5{ri - #*•) 



df/\(Ar,f,f/)[T(fr) - T(f)]p (2 \f, ft) 



(3.20) 
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Expanding now the r/-dcpcndcncc around r in powers of Ar and collecting the first significant order, the 



right hand side of eq. (3.20) becomes 



V ^ {~^~ Jo ^ Ar ' 4 '^ Ar ^^ 9 ^ Ar '^ f ^] VT (3.21) 



Thus, the macroscopic thermal conductivity takes the form 



£ = /9 2 |y^ dArAr 4 X(Ar,f,r)g(Ar,r,f)^=p 2 X A (3.22) 



where the coefficient A4 is analogous to 774 defined in eq. ( 3.11 ). This result has been given in eq.(3.25) 
ofref. @. 



The total energy transport equation is obtained by adding eqs. (3.14) and (3.18), yielding 



J^pe = -V- ^mpv + pe) - V-pv + V- (v-a^j + V • £VT (3.23) 

where e = e m + e u and we have introduced the macroscopic internal energy density of the DPDE particles 
system, e, as the sum of the internal kinetic energy, the interaction potential energy between the particles 
and the internal energy of the particles 

pe = ^kTp + p«f + pe u (3.24) 



which is the same equation of state as given in eq. (3.12) of ref. 0. Note that Eq. ( 3.23| ) corresponds to 



a transport equation for a conserved magnitude. This form reflects the fact that the total energy of the 
system is conserved in the particle-particle interaction. Finally, making use of eq. ( |3.8| ) to eliminate the 
macroscopic kinetic energy transport, one arrives at the transport equation for the macroscopic internal 
energy 

d 

—pe = -V • pev-pV -v + a : W+ V • £VT (3.25) 
at 

Each term on the right hand side of this equation has an interpretation in the framework of the macro- 
scopic transport phenomena. The first term stands for the advection of the macroscopic internal energy 
due to macroscopic fluid motion. The second and third terms are, respectively, a change in the internal 
energy due to the work due to an expansion of the fluid, and the viscous heating of the system due to 
viscous forces. From the fourth and last term we can define the macroscopic heat flow 

J q = -£VT (3.26) 

which is precisely the macroscopic expression of the Fourier law. 

Although the theoretical approach developed here is only applicable to the large friction limit, we end 
this section by including a qualitative evaluation of the kinetic contributions to the thermal conductivity, 
to enable a comparison with the simulation data. 

Effectively, we can estimate the time t p taken by a particle to loose the memory of its initial momentum 
as being given by 

(3 - 27) 



as follows from a balance between dpi/dt ~ Pi/t p and the force given in eq. (2.1), times the number of 
particles inside the interaction range, pr^. Since the particles have an average velocity v ~ yJkT '/m, we 
can estimate the momentum penetration depth as 
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On the other hand, we can also estimate the characteristic relaxation time for the internal energy of a 
given particle (much larger than the kinetic energy since we assume that <fi 3> k) from a balance between 



dui/dt ~ 4>T/t q and the mesoscopic heat flow, according to eq. (2.3) This leads to 



t q ~ (3.29) 

where we have also considered that pr\ particles are simulatenously interacting. Hence, if t q <C t p the 
particle covers a distance l q given by 

in the time the energy is dissipated. Thus, if l q <C l p , the energy transport is dominated by the dissipative 
contribution and the energy is delivered before the momentum has relaxed. Thus, the contribution to 
the energy flow due to the motion of the particles is given by p6| 

J q ~ cj)pvl q X7T (3.31) 

Therefore, in the region dominated by the dissipative interactions, the kinetic contribution to the thermal 
conductivity is subdominant and has the functional form 

C K i ~ p4>vl q ~ (3.32) 

Notice that Cki scales as 1/Lo, according to our initial assumption. A similar reasoning leads to the 
1 /(^-dependence of the kinetic contribution to the viscosity in the region where the momentum tr anspor t 
is dominated by the friction between particles, as has been already obtained by different authors fl3| . |27|| . 

If Lq is reduced, we will eventually have that l q > l p . Clearly, in this region the energy transport 
is dominated by the kinetic contribution, which is much more effective than the direct transmission of 
heat between the particles. It is crucial to realize, however, that the flow of particles is dominated by its 
diffusive motion instead of the inertial motion of the previous domain, since the momentum has had time 



to relax. The characteristic relaxation time is still given by eq. (3.29), since the mechanism for energy 
delivery keeps on being the direct heat transport between particles. However, l q is calculated from the 
mean displacement of a brownian particle, that is 



U^S (3.33) 

V Co/Wf 

where the diffusion coefficient has been estimated from a Stokes- Einstein relation D ~ kT/i^opr^. Thus, 
the characteristic velocity of this diffusive flow is l q /t q , from which we obtain the energy flow in this 
regime 

J q ~ (j>p l -?-l q VT ~ (t> P DVT (3.34) 

tq 

Hence, from this expression we infer the kinetic contribution to thermal conductivity which, in addition, 
will be the leading contribution in this region 
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Ckii ~ 4>pD ~ (3.35) 
Cor? 



The preceeding reasoning leads to a resulting thermal conductivity independent of the parameter Lq in 
the region where the dissipative contribution to the thermal conductivity is small. In this respect, this 
new result is essentially different from that reported for the viscosity coefficient in the isothermal DPD 
model lilJll. 



IV. SIMULATION RESULTS 



The particular system used to investigate the transport properties of the DPDE model consists of 
N = 10000 particles in a cubic box of periodic boundary conditions in all three dimensions, so that no 
walls are considered. The lateral size of the box has been set to A 1 / 3 in dimensionless units, and hence 
the number density p* is always unity. We have studied essentially dense systems, in the sense that 
many particles interact with each other at one time, and therefore, the interaction ranges are larger than 
l/p 1 / 3 . The choice of p* = 1 also implies that the momentum penetration depth is smaller than, or of 
the order of, the mean distance between particles. 

Before taking any measurment, the system is allowed to reach thermal equilibrium. According to the 
expressions given in eqs. (3.11) and ( |3.22 ), and the form chosen in eqs. (2.9) and (2.10) for the dissipative 



coefficients, the predicted values of the thermal conductivity and the shear viscosity are, respectively 

£ = lf^ <«» 

l-^*b (4-2) 

Notice the resulting temperature dependence on the thermal conductivity, which arises from the depen- 
dence of the mesoscopic thermal conductivity on the particles' temperatures. 

The measurement of the transport properties in the system has been carried out by tracking the 
relaxation of perturbations externally induced in the system. The hydrodynamic behavior predicted for 
the DPDE system in the prec ee ding secti on su ggests that small perturbations should relax according to 
linearised versions of eqs. (|J), |[|) and ( EH )- In particular, one finds that a given Fourier component 



of wavevector k of a temperature perturbation field in the system should relax as 

6T k (t) = 5T k (0)e- Kk ' 2t (4.3) 



where k = C/c v p is the thermal diffusivity, with 

where c v is the constant volume heat capacity per particle (related to the value given in eq. (3.22) of 
ref. H according to c v = C v /N). Notice that c v ^ <j), due to the contribution of the energy stored in 
other degrees of freedom different from the internal energy of the particle. Superposed to this relaxation 
there is also the effect of the temperature fluctuations induced by the random forces and random heats. 
However, for small values of B, the amplitude of the random temperature fluctuations is very small and, 
thus, can be ignored. 

We have obtained the relaxation of temperature fluctuations, c v p5T k , as the difference between per- 
turbations in the internal energy density S(pe)k, externally induced in the system at t = 0, and c v T5pk, 
which are the directly observable variables. To obtain the value of the thermal conductivity, we have 
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fitted an exponential decay to the curves obtained in each simulation. For all the simulations dominated 
by the dissipative contribution, we get clear exponential decays, with regression coefficients typically 
of 0.999 or higher, and error bars for the measured exponent of about 1% of its value. In the region 
dominated by the kinetic contribution, slight deviations from a purely exponential decay are observed. 
However, correlation coefficients typically of 0.98 and 0.99 are still found, and the error is about 10%. 
The measured values for the thermal conductivity as given by means of the procedure described above 
have been compared with those obtained from the application of a temperature gradient to the system 
and measuring the resulting heat flow through the system once the steady state was reached. The results 
obtained by both methods are in agreement. 

We have performed two series of measuremets of the thermal conductivity for two different radii r\. 
The interaction number is then given by 

An , An . ... , 
n = ypri = y Af (4-5) 

We have plotted the results for C* = Cml/(j)(o against the dimensionless parameter X = Cr* x p* 2 /T* 2 , 
for values of X ranging from 1 to 1000, thus covering three decades of this parameter. The results are 



then compared with the theoretical prediction, given in eq. (4.1), that reads 



in dimensionless form. The dimensionless temperature was about 1.1 and the density was unity. The 
range of the interaction r* x was set to 1.24 and 2.88, so that the interaction number was n = 8 and 
n = 100, respectively. The dimensionless shear viscosity coefficient takes the form 

* ' 27r 5 2 

11 = % = i^ rc p (47) 

According to our set of dimensionless variables, the viscosity coefficient depends only on the density and 
interaction range and, hence, is constant for each set of simulations characterized by a given value of n. 
We have measured a viscosity coefficient if = 0.363 for the case n = 8 and 77* = 0.894 for n = 100. We 
set the parameter B to 10 -5 to avoid negative values of the energy during the simulation for large C 
values. Likewise, the time-step was chosen to be St* = 10~ 2 or 5t* = 10~ 3 in the critical cases regarding 
the possibility of negative values of the particles' energy. 

In fig. 1 we plot our results for the thermal conductivity as a function of the parameter X. For 
both values of n we observe a linear dependence of C* in this parameter, for values of X larger than 
315/2-7T ~ 50, which roughly corresponds to the point at which the thermal conductivity is about 1. Thus, 
we can conclude that the dissipative contribution dominates the thermal conductivity in this region, and 



that the dependence in X is well captured by eq. (4.6). In addition, for the case n = 100, we find a much 



better agreement between the theoretical predictions and the simulation results than in the case n = 8, 
whose value is roughly half of the theoretical one in all of the range. This is a rather unexpected result 
which is caused neither by the particle number fluctuations in the sample nor by significant deviations 
of g(r) from 1, as we have checked. Such a behavior is equally independent of the way of measuring 
the property, since good agreement has been found between the data obtained from both the relaxation 
method and heat flow in a temperature gradient. Furthermore, data reported in ref. |l7| for a system 
of frozen DPDE particles also show the same kind of dependence of the thermal conductivity with n, so 
that it cannot be attributed to a possible coupling between energy and momentum transport. This effect 
has no equivalent in the case of the viscosity coefficient which, in our dimensionless variables, varies only 
with the range of the interaction at constant p* . 

In fig. 2 we focus our attention on the range < X < 50, so that the dissipative contribution to the 
thermal conductivity is subdominant. We plot the results for the two series of data (n = 8 and n = 100) 
as a function of the parameter Y = T* 3 /Cr^ 3 ~ £-ki, to compare with the prediction given in eq. fl3.32p . 
Both series of data agree with a thermal conductivity independent of Lq (C in dimensionless variables) 
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in this range of values of Lq. Althoug h the data are less reliable than in the previous range, we can 



consider that the behavior given in eq. ( 3.35 ) is confirmed and that the transport of energy for small Lq 
is due to particle diffusion. We should point out, however, that the data with n — 8 show an increase of 
the thermal conductivity in passing from the dissipative regime to the diffusive regime, while the data 
for n — 100 shows a smooth crossover. This can attributed to the fact that for n = 8 the crossover 



between both regimes passes first through the behavior described in eq. ( 3.32 ). In the n = 100 case, the 
dissipative and diffusive regimes simply overlap in the crossover. 

To end this section, we report data on the shear viscosity coefficient obtained for the model. According 
to our choice of parameters, the shear viscosity is only a function of the range of interaction r^. Thus, we 
have performed simulations and plotted in fig. 3 the viscosity coefficient as a function of the parameter 
Z = p* 2 rt 5 . As expected, we obtain a good agreement between the theoretical prediction given in eq. 



(4.7) for values of the parameter Z larger than 1575/2-7T ~ 250, which roughly corresponds to a viscosity 
rj* ~ 1. In this region, however, a slight curvature can be observed. For values of Z smaller than 250, as 
is to be expected, the momentum transport is dominated by kinetic effects. Here, we observe the same 
kind of qualitative behavior as found for the isothermal DPD model as has been discussed elsewhere 
~!l6|]i~5|. 



V. CONCLUSIONS 



In this paper we have analyzed the transport properties of the new DPDE model from both, a simple 
theoretical approach and computer simulations. This model is addressed at the simulation of fluctuating 
fluids in which momentum as well as heat transport are involved, and has already shown well defined 
equilibrium properties ||. In this respect, the model is complete in the sense that it is thermodynamically 
consistent and the five hydrodynamic fields can be correctly described. The complete understanding of 
its features and the refinement of its capabilities, however, still demands a great deal of additional effort. 

After a brief review of the algorithm, we have first of all derived the dynamic properties of the DPDE 
model based on a local equilibrium assumption. In this framework, we have obtained the transport 
equations and approximate expressions for the transport coefficients which here can be functions of the 
temperature. In this approach, it is implicitly assumed that an initial probability distribution for the 
relevant variables of the system will relax to the local equilibrium form after a time t p for the momentum 
and t q for the energy, according to eqs. ( 3.27j ) and ( 3.29| ), respectively. These characteristic times must 



be much smaller than those related to the changes in the hydrodynamic fields for the whole approach 
to be valid. We have qualitatively argued that the local equilibrium is the dominant contribution in an 
expansion in inverse powers of the dissipative coefficients Lq and Co, so that the transport coefficients 
can be expressed as series expansions of this particular form. This is also the case in the derivation of 
the Smoluchowski equation from the Fokker-Planck equation |i~8|-|20[]. Therefore, in the range of validity 
of this point of view, the kinetic contributions are always subdominant. The analytic form for these 
subdominant kinetic contributions has been reported in the literature (l^j2^] in the case of the viscosity, 
and are given in eq. ( HH ) in the case of the thermal conductivity. 



Second, we have also studied the behavior of the system when the dissipative contribution becomes 
comparable or even smaller than the kinetic contribution to the dissipative coefficients. In this regime, the 
theoretical picture based on the local equilibrium assumption fails. It is important to realize, however, 
that the form of the series expansion of the transport coefficients in inverse powers of Lq and Co will 
change towards another functional dependence. This is in fact one of the more important findings of this 
paper which also applies in the case of viscosity in the isothermal DPD model. In the case of the thermal 
conductivity we have found that t he th ermal conductivity at small values of Lq should be driven by 



particle diffusion, according to eq. ( 3 . 3E ) . We find that C is independent of Lq and that the simulations 
performed seem to agree with this interpretation. Thus, the increase of the thermal conductivity at small 
values of Lq is attributed to a crossover from the regime dominated by the dissipative interactions to a 
regime dominated by diffusive transport. This crossover region depends on the interaction number n. 

Third, we have shown the first simulation results for both, the shear viscosity coefficient and the 
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thermal conductivity for the DPDE model. Qualitatively, the shear viscosity for the DPDE behaves in 
the same way as that of the isothermal DPD model, as could have been expected. Furthermore, we have 
explored different parameter regions that cover only part of all the possibilities of the DPDE model. 
Particular attention has been paid to non-dilute systems, aiming at a reliable comparison between the 
simulation data and the theoretical results in the large friction limit. We have also reported results in 
regions where the theoretical assumptions fail and new qualitative behaviors have been found for the 
thermal conductivity. 

Fourth and last, an important result of our simulation analysis is the difference (of about a factor of 
two between the two series of data) observed in the thermal conductivity, in the region of large Lq, for 
the two values of the interaction number analyzed. We find no qualitative explanation for this difference, 
which becomes larger as the number of particles inside the interaction range diminishes. The origin of this 
discrepancy could be the same as that found for the viscosity in the isothermal DPD model, according 
to ref. The fact that the CPU time grows with the interaction number, makes the regime of small 
n be of great interest by itself. 
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FIGURE CAPTIONS 



Fig. 1: Dimensionless thermal conductivity as a function of X = Cr*^ p* 2 /T* 2 . Squares stand for data 
for n = 8 while circles represent data for n = 100. The solid line is the theoretical prediction given in 
eq.Q 

Fig. 2: Dimensionless thermal conductivity as a function of Y = T* 3 /Cr^ 3 . Squares stand for data for 
n = 8 while circles represent data for n = 100. The inset shows a detailed view of the n = 8 data. 

Fig. 3: Dimensionless shear viscosity as a function of Z = p* 2 rt 5 . The solid line stands for the theoretical 



prediction given in eq. (4.7), while the circles indicate simulation data points. 
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